Content under a Creative Commons Attribution license, CC-BY 4.0; code under MIT license. (c)2014 Lorena A. Barba, Olivier Mesnard. Thanks: NSF for support via CAREER award #1149784.
Version 0.3 -- March 2014

Method of images

In the previous AeroPython lessons, you've seen how the fundamental solutions of potential flow can be combined to produce stream-line patterns that represent flow around objects—e.g., a Rankine oval and a circular cylinder.

You may ask now: can we get flow around a straight wall? Yes, we can!

The method of images is simply this: put a singularity near a "wall" by adding an "image" singularity on the other side of the wall. Sometimes, this is also referrered to as aerodynamic interferences.

In addition to that, this notebook also gives us the opportunity to introduce the notion of classes in Python. This is a very useful way to organize your code, that becomes crucial as programs get more complex.

As usual, we will start by importing the libraries and creating the computational grid. Let's get that out of the way.


In [1]:
import numpy
import math
from matplotlib import pyplot

In [2]:
N = 50                                # Number of points in each direction
x_start, x_end = -2.0, 2.0            # x-direction boundaries
y_start, y_end = -1.0, 1.0            # y-direction boundaries
x = numpy.linspace(x_start, x_end, N)    # computes a 1D-array for x
y = numpy.linspace(y_start, y_end, N)    # computes a 1D-array for y
X, Y = numpy.meshgrid(x, y)              # generates a mesh grid

Source near a plane wall

If a source is located near a wall, the flow pattern will be distorted by the wall. Imagine a source located at $y=y_\text{source}$ next to a wall at $y=0$. The boundary condition at the wall requires the flow to be tangent there—for a horizontal wall, $v=0$. The effect of the wall on the source flow is identical to placing another source (an image) of equal strength, located at $y=-y_\text{source}$.

Ok, now it is time to play smart and avoid repeating code that is wasting our precious time! In a previous notebook, we already introduced the notion of function in Python, and now we are going to move forward with the creation of classes.

A class is a bundle of data (parameters and variables) and "methods" or functions that work with this data. It's a very tidy way of organizing code. As we create longer and more complex codes, this helps us manage complexity. The code becomes easier to maintain, modify and extend.

We will define a class called Source that will contain information related to a source. A singularity such as a source is defined by its strength and its location in the domain. Therefore, our class Source will have three attributes defining it:

  • strength: the strength of the source.
  • x: the location of the source on the $x$-axis.
  • y: the location of the source on the $y$-axis.

What would we like to do after the definition of our source? We want to compute the velocity field, as well as the stream function. So, in our class Source, we have to implement two methods (functions included in a class): one to compute the velocity (a function called velocity) and the other one to calculate the stream function (a function called stream_function).

Let's have a detailed look of the methods:

  • method velocity: is a function used to calculate the velocity on the mesh (X,Y) due to a source. Therefore, our method will have two arguments, X and Y, and generate two new attributes of our source: u and v, the velocity components of the source.
  • method stream_function: is a function used to calculate the stream function on the mesh (X,Y) due to a source. Therefore, our method will have two arguments, X and Y, and generate one new attribute of our source: psi, the stream function of the source.

In addition to these two, every class needs to have a constructor, which is simply a fancy word for initializing the data. It is a function that is always called __init__:


In [3]:
class Source:
    """Contains information related to a source (or sink)."""
    def __init__(self, strength, x, y):
        """Initializes the singularity.
        
        Arguments
        ---------
        strength -- strength of the singularity.
        x, y -- coordinates of the singularity.
        """
        self.strength = strength
        self.x, self.y = x, y
    
    def velocity(self, X, Y):
        """Computes the velocity field generated by the singularity.
        
        Arguments
        ---------
        X, Y -- mesh grid.
        """
        self.u = self.strength/(2*math.pi)*(X-self.x)/((X-self.x)**2+(Y-self.y)**2)
        self.v = self.strength/(2*math.pi)*(Y-self.y)/((X-self.x)**2+(Y-self.y)**2)
    
    def stream_function(self, X, Y):
        """Computes the stream-function generated by the singularity.
        
        Arguments
        ---------
        X, Y -- mesh grid.
        """
        self.psi = self.strength/(2*math.pi)*numpy.arctan2((Y-self.y), (X-self.x))

What's this self business? When we call the class Source, it will look like a function call. For example: source = Source(1, 0, 0), to create a source of strength 1 located at the origin. Python automatically calls the constructor function with the parameters (self, 1, 0, 0), i.e., Python adds self to the list of parameters, to mean "the source that needs to be created now with the parameters that follow."

Once we've created a source, we can call its two methods to compute the velocity field and the stream function on the mesh (X,Y). Like when using libraries, we'll indicate the method of a class by using the dot notation (as shown below).

The Cartesian velocity components are given by:

$$u = \frac{\sigma}{2\pi}\frac{x-x_\text{source}}{\left(x-x_\text{source}\right)^2+\left(y-y_\text{source}\right)^2}$$$$v = \frac{\sigma}{2\pi}\frac{y-y_\text{source}}{\left(x-x_\text{source}\right)^2+\left(y-y_\text{source}\right)^2}$$

and the stream function is defined by:

$$\psi = \frac{\sigma}{2\pi}\arctan\left(\frac{y-y_\text{source}}{x-x_\text{source}}\right)$$

Now, let's see how this works.


In [4]:
strength_source = 1.0                  # strength of the source         
x_source, y_source = 0.0, 0.5          # location of the source

# creates a source (object of class Source) 
source = Source(strength_source, x_source, y_source)

# computes the velocity field and the stream-function on the mesh grid
source.velocity(X, Y)
source.stream_function(X, Y)

Notice the dot notation? With the dot, we are saying: "take the velocity function of the object source that we just created, and execute it." You see, the velocity function is associated to the data of the particular source, because it is part of a class.

The image of the source will also be created using the class Source, with the same strength but different location. You probably start to see the efficiency of having created a class.


In [5]:
# creates the image of the source and computes velocity and stream-function
source_image = Source(strength_source, x_source, -y_source)
source_image.velocity(X, Y)
source_image.stream_function(X, Y)

Using superposition, we can compute the streamlines of a source in the vicinity of a wall. The superimposition of the two sources leads to the following velocity field:

$$u = \frac{\sigma}{2\pi} \left( \frac{x-x_\text{source}}{\left(x-x_\text{source}\right)^2+\left(y-y_\text{source}\right)^2} + \frac{x-x_\text{source}}{\left(x-x_\text{source}\right)^2+\left(y+y_\text{source}\right)^2} \right)$$$$v = \frac{\sigma}{2\pi} \left( \frac{y-y_\text{source}}{\left(x-x_\text{source}\right)^2+\left(y-y_\text{source}\right)^2} + \frac{y+y_\text{source}}{\left(x-x_\text{source}\right)^2+\left(y+y_\text{source}\right)^2} \right)$$

and the following stream function:

$$\psi = \frac{\sigma}{2\pi} \left( \arctan\left(\frac{y-y_\text{source}}{x-x_\text{source}}\right) + \arctan\left(\frac{y+y_\text{source}}{x-x_\text{source}}\right) \right)$$

In [6]:
# superposition of the source and its image
u = source.u + source_image.u
v = source.v + source_image.v
psi = source.psi + source_image.psi

# plots the streamlines
%matplotlib inline

size = 10
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(source.x, source.y, color='#CD2305', s=80, marker='o')
pyplot.scatter(source_image.x, source_image.y, color='#CD2305', s=80, marker='D')
pyplot.axhline(0., color='k', linestyle='--', linewidth=4);


Vortex near a plane wall

In a similar manner, we write another class called Vortex, requiring as data its strength and location. The class will have two methods: one to compute the velocity of the vortex and the other one to calculate the stream-function.

But we also always need a constructor, a function called __init__, to initialize the data of a vortex:


In [7]:
class Vortex:
    """Contains information related to a vortex."""
    def __init__(self, strength, x, y):
        """Initializes the vortex.
        
        Arguments
        ---------
        strength -- strength of the vortex.
        x, y -- coordinates of the vortex.
        """
        self.strength = strength
        self.x, self.y = x, y
        
    def velocity(self, X, Y):
        """Computes the velocity field generated by a vortex.
        
        Arguments
        ---------
        X, Y -- mesh grid.
        """
        self.u = +self.strength/(2*math.pi)*(Y-self.y)/((X-self.x)**2+(Y-self.y)**2)
        self.v = -self.strength/(2*math.pi)*(X-self.x)/((X-self.x)**2+(Y-self.y)**2)
        
    def stream_function(self, X, Y):
        """Computes the stream-function generated by a vortex.
        
        Arguments
        ---------
        X, Y -- mesh grid.
        """
        self.psi = -self.strength/(4*math.pi)*numpy.log((X-self.x)**2+(Y-self.y)**2)

It is now very quick to create two object of type Vortex for a new vortex and its image.


In [8]:
strength_vortex = 1.0                  # strength of the vortex
x_vortex, y_vortex = 0.0, 0.5          # location of the vortex

# creates a vortex and computes velocity and stream-function
vortex = Vortex(strength_vortex, x_vortex, y_vortex)
vortex.velocity(X, Y)
vortex.stream_function(X, Y)

# creates the image of the vortex and computes velocity and stream-function
vortex_image = Vortex(-strength_vortex, x_vortex, -y_vortex)
vortex_image.velocity(X, Y)
vortex_image.stream_function(X, Y)

Applying the principle of superposition, we are able to get the streamlines representing a vortex near a wall.


In [9]:
# superposition of the vortex and its image
u = vortex.u + vortex_image.u
v = vortex.v + vortex_image.v
psi = vortex.psi + vortex_image.psi

# plots the streamlines
size = 10
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(vortex.x, vortex.y, color='#CD2305', s=80, marker='o')
pyplot.scatter(vortex_image.x, vortex_image.y, color='#CD2305', s=80, marker='D')
pyplot.axhline(0., color='k', linestyle='--', linewidth=4);


Motion of a vortex pair near the ground

We're having too much fun. Let's get more singularities near a wall. Two vortices will make you dizzy ...


In [10]:
strength_vortex = 1.0                  # absolute value of each vortex strength
x_vortex1, y_vortex1 = -0.1, 0.5       # location of the first vortex
x_vortex2, y_vortex2 = +0.1, 0.5       # location of the second vortex

# creates two vortices at different locations
vortex1 = Vortex(+strength_vortex, x_vortex1, y_vortex1)
vortex2 = Vortex(-strength_vortex, x_vortex2, y_vortex2)

# computes the velocity and stream-function for each vortex
vortex1.velocity(X, Y)
vortex1.stream_function(X, Y)
vortex2.velocity(X, Y)
vortex2.stream_function(X, Y)

# creates an image for each vortex
vortex1_image = Vortex(-strength_vortex, x_vortex1, -y_vortex1)
vortex2_image = Vortex(+strength_vortex, x_vortex2, -y_vortex2)

# computes the velcoity and stream-function of each image
vortex1_image.velocity(X, Y)
vortex1_image.stream_function(X, Y)
vortex2_image.velocity(X, Y)
vortex2_image.stream_function(X, Y)

In [11]:
# superposition of the vortex pair and its image
u = vortex1.u + vortex2.u + vortex1_image.u + vortex2_image.u
v = vortex1.v + vortex2.v + vortex1_image.v + vortex2_image.v
psi = vortex1.psi + vortex2.psi + vortex1_image.psi + vortex2_image.psi

# plot the streamlines
size = 10
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(vortex1.x, vortex1.y, color='#CD2305', s=80, marker='o')
pyplot.scatter(vortex2.x, vortex2.y, color='g', s=80, marker='o')
pyplot.scatter(vortex1_image.x, vortex1_image.y, color='#CD2305', s=80, marker='D')
pyplot.scatter(vortex2_image.x, vortex2_image.y, color='g', s=80, marker='D')
pyplot.axhline(0., color='k', linestyle='--', linewidth=4);


Doublet near a plane wall parallel to a uniform flow

Last one ... a doublet next to a wall. We need a new class, this time named Doublet. Same as before!


In [12]:
u_inf = 1.0    # free-stream speed

# calculates the velocity and stream-function of the free-stream flow
u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)
psi_freestream = u_inf * Y

In [13]:
class Doublet:
    """Contains information related to a doublet."""
    def __init__(self, strength, x, y):
        """Initializes the doublet.
        
        Arguments
        ---------
        strength -- strength of the doublet.
        x, y -- coordinates of the doublet.
        """
        self.strength = strength
        self.x, self.y = x, y
        
    def velocity(self, X, Y):
        """Computes the velocity field generated by a doublet.
        
        Arguments
        ---------
        X, Y -- mesh grid.
        """
        self.u = -self.strength/(2*math.pi)*\
                ((X-self.x)**2-(Y-self.y)**2)/((X-self.x)**2+(Y-self.y)**2)**2
        self.v = -self.strength/(2*math.pi)*\
                2*(X-self.x)*(Y-self.y)/((X-self.x)**2+(Y-self.y)**2)**2
            
    def stream_function(self, X, Y):
        """Computes the stream-function generated by a doublet.
        
        Arguments
        ---------
        X, Y -- mesh grid.
        """
        self.psi = -self.strength/(2*math.pi)*(Y-self.y)/((X-self.x)**2+(Y-self.y)**2)

In [14]:
strength_doublet = 1.0                # strength of the doublet
x_doublet, y_doublet = 0.0, 0.3       # location of the doublet

# creates a doublet (object of class Doublet)
doublet = Doublet(strength_doublet, x_doublet, y_doublet)

# computes the velocity and stream-function of the doublet on the mesh
doublet.velocity(X, Y)
doublet.stream_function(X, Y)

# creates the image of the doublet
doublet_image = Doublet(strength_doublet, x_doublet, -y_doublet)

# computes the velocity and stream-function of the image on the mesh
doublet_image.velocity(X, Y)
doublet_image.stream_function(X, Y)

In [15]:
# superposition of the doublet and its image to the uniform flow
u = u_freestream + doublet.u + doublet_image.u
v = v_freestream + doublet.v + doublet_image.v
psi = psi_freestream + doublet.psi + doublet_image.psi

# plots the streamlines
size = 10
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(doublet.x, doublet.y, color='r', s=80, marker='o')
pyplot.scatter(doublet_image.x, doublet_image.y, color='r', s=80, marker='D')
pyplot.axhline(0., color='k', linestyle='--', linewidth=4);



Please ignore the cell below. It just loads our style for the notebook.

In [16]:
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom.css', 'r').read()
    return HTML(styles)
css_styling()


Out[16]: